% Edited by Natasha Savage 30Nov10

% cargo "sinks" into cells destined for endocytosis

    if ( endo_cell.num > 0 )
      %  endo_time includes that 15-s endo_red_time
      %  find "green" bins
      ind = find(endo_cell.endo_time < 0);        
     
     [ind2,Xi,Yi,interpD] = sink_ind_02g (endo_cell.posx(ind),endo_cell.posy(ind), dx, N);      
     
      [endo_cell.vSNARE(ind), vSNARE] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.vSNARE(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),vSNARE,dx,endo_cell.dx(ind),dt,Dconst,1);
      [endo_cell.Rec(ind), Rec] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.Rec(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),Rec,dx,endo_cell.dx(ind),dt,Dconst,0); %Receptor alone does not sink %Masha dev dec2012
      [endo_cell.RecA(ind), RecA] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.RecA(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),RecA,dx,endo_cell.dx(ind),dt,Dconst,1); %Receptor bound to pheromone sinking %Masha dev oct2012
      
      
      [endo_cell.Cdc42T(ind), Cdc42T] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.Cdc42T(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),Cdc42T,dx,endo_cell.dx(ind),dt,Dconst,0);
      [endo_cell.Cdc42D(ind), Cdc42D] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.Cdc42D(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),Cdc42D,dx,endo_cell.dx(ind),dt,Dconst,0);
      [endo_cell.BemGEF42(ind),  BemGEF42] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.BemGEF42(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),BemGEF42,dx,endo_cell.dx(ind),dt,Dconst,0);
      [endo_cell.GDI42(ind), GDI42] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.GDI42(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),GDI42,dx,endo_cell.dx(ind),dt,Dconst,0); 
      
      [endo_cell.BemGEF(ind), BemGEF] = ...
          sink_02g (ind2,Xi,Yi,interpD, ...
          endo_cell.BemGEF(ind),endo_cell.posx(ind),...
          endo_cell.posy(ind),BemGEF,dx,endo_cell.dx(ind),dt,Dconst,0);                   
       
    end 
ticd=tic;
    % diffuse
    % membrane
    vSNARE   = reshape(Hop\reshape(vSNARE,  N^2,1),N,N);
    Rec      = reshape(Hop\reshape(Rec,     N^2,1),N,N); %Added for The Receptor Complex %Masha dev oct2012
    RecA     = reshape(Hop\reshape(RecA,    N^2,1),N,N); %Added for The Receptor Complex %Masha dev oct2012
    Cdc42T   = reshape(Hop\reshape(Cdc42T,  N^2,1),N,N);
    BemGEF42 = reshape(Hop\reshape(BemGEF42,N^2,1),N,N);
    BemGEF   = reshape(Hop\reshape(BemGEF,  N^2,1),N,N);
    Cdc42D   = reshape(Hop\reshape(Cdc42D,  N^2,1),N,N);
    GDI42    = reshape(Hop\reshape(GDI42,   N^2,1),N,N);
    
    % calculate new cytoplasmic concns by averaging over all cells
    BemGEFc(:,:) = mean(mean(BemGEFc));
    GDI42c(:,:)  = mean(mean(GDI42c));
    GDIc(:,:)    = mean(mean(GDIc));    
tdiffuse=[tdiffuse toc(ticd)];
 %% exocytosis
 
 % Cable attach
 
 if length(cable)/2 < MAXcable                           % attach new cables
     cable_add = [];                                    % store new cable positions
      for  i = 1:(MAXcable - (length(cable)/2))
          if random('unif',0,1) <= attachPROB           % attach a cable

              dice = random('unif',0,1,N,N).*ves_prob;  % random distribution within TOP bins
              [indx,indy] = find(dice==max(max(dice))); % [indx,indy] = index of winning bin, bottom-left of 2x2 spot

              if (length(indx)>1) % in case 2 bins get same random number          
                indx = indx(1)
                indy = indy(1)
              end
              % to choose the corner, bottom left default
              % top or bottom
              if rand()>0.5    % choose top
                indy = indy-1;
                if indy == 0; indy = N; end
              end
              % left or right
              if rand()>0.5    % choose right
                indx = indx-1;
                if indx == 0; indx = N; end
              end
              
              cable_add = [cable_add, indx, indy];      % store new cable positions
          end
      end
 end
      
% Cable detatch
if length(cable) ~= 0                                  % detatch cables
    cable_rmv = [];
    for  i = 1:2:length(cable)
        if random('unif',0,1) <= detatchPROB            % detatch a cable
            cable_rmv = [cable_rmv,i,i+1];
        end
    end
    cable(cable_rmv) = [];
end

cable = [cable cable_add];
cable_add = [];

% cs2 = [NaN(1,N)];
% cs2(1:length(cable)) = cable;

    %% an exocytotic event may happen

   dice = random('unif',0,1); % dice=random number between 0 and 1  
   if ( dice < dt/periodIn && curr_t && length(cable)>0) % exocytosis will happen if dice<dt/periodIn (=0.05/2.4=0.02: 2% chance)  
      exo_cnt = exo_cnt + 1;

      % Pick cable
      dice = random('unif',0,1,1,(length(cable)/2));        
      i    = find(dice==max(dice));
      if (mod(i,2)==0) %if the index is 2th in a found pair (Y coordinate of a cable), shift to the X index.
          i=i-1;
      end
      %i = (2*i)-1;
      indx = cable(i);
      indy = cable(i+1);  

      % Radially 'stretch' [Cdc42] distribution, set up variables for interpolation

      r_exo  = nVInlen*dx/(pi^0.5);  % radius of area added - old membrane size
                                     % area added is (nVInlen*dx)^2      
      C = N*dx;
      hole_radial_exo_area_05;                      % call script 'hole' to compute

      % transform all concentrations so vesicle event is in the centre
      vesicle_position_center;
      
      % add naked membrane to external membrane  
      vSNARE   = exo_stretch_radial_01 ( vSNARE,  Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2)); 
      Rec = exo_stretch_radial_01 ( Rec,Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2)); %Added for Receptor Complex %Masha dev oct2012
      RecA     = exo_stretch_radial_01 ( RecA,    Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2)); %Added for Receptor bound to pheromone %Masha dev oct2012
      Cdc42T   = exo_stretch_radial_01 ( Cdc42T,  Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2));
      BemGEF42 = exo_stretch_radial_01 ( BemGEF42,Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2));
      BemGEF   = exo_stretch_radial_01 ( BemGEF,  Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2));
      Cdc42D   = exo_stretch_radial_01 ( Cdc42D,  Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2));
      GDI42    = exo_stretch_radial_01 ( GDI42,   Fwarp,  x0, y0, x1, y1,[N/2,(N/2)+1],[N/2,(N/2)+1],N^2/(N^2+nVInlen^2));
      clear xold yold        
            
      % transform all concentrations back
      vesicle_position_restore;
      
      %% move endocytic patches after exocytosis

      dx_new  = dx*sqrt(N^2+nVInlen^2)/N;
      C = N*dx_new;                  % move patches on new membrane 
      endo_cell = moving_patches_exo_radial(r_exo,(indx*dx_new)+(dx_new/2),(indy*dx_new)+(dx_new/2),endo_cell,C);                               
 %%
       
      % adjust size of external membrane, and diffusion operator
      dx  = dx*sqrt(N^2+nVInlen^2)/N;
      eta = ((Rnew_mult^3 * dx^3) - (((Rnew_mult * dx)-mem_depth)^3))/cyt_mult;
      Dxx = Dxx*(N^2+nVInlen^2)/N^2;
      Hop = speye(N^2) - kron(speye(N),Dxx) - kron(Dxx,speye(N));                

      % adjust size of internal membrane
      dx2new = dx2*sqrt(N^2-(nVInlen*dx/dx2)^2)/N;
      
      if (indx<=0 | indx > N | indy <=0 | indy>N)
        display(['dice= ' num2str(dice) 'indx = ' num2str(indx) ' indy = ' num2str(indy)]);
      end

      
      % Set [cargo] in new patch of membrane 
      vSNARE([indx,mod(indx,N)+1],[indy,mod(indy,N)+1]) = 10*vSNAREic; %Masha dev oct2012 %pbc wrap
      GDI42( [indx,mod(indx,N)+1],[indy,mod(indy,N)+1]) = GDI42ic;
      Cdc42D([indx,mod(indx,N)+1],[indy,mod(indy,N)+1]) = Cdc42Dic;
      Rec([indx,mod(indx,N)+1],[indy,mod(indy,N)+1]) = 10*Recic; % Bulk traffic of the Receptor %Masha dev oct2012
      
%      fprintf(rec_count_time,'%5.4f %5.4f ',curr_t, 10*Recic);
     
      %remove transferred cargo from internal compartment
      vSNAREic = vSNAREic*(N^2*dx2^2-40*dx^2)/(N^2*dx2new^2);  
      Recic= Recic*(N^2*dx2^2-40*dx^2)/(N^2*dx2new^2); %Masha dev sept2012
      dx2 = dx2new;    
      
      exotimes =[exotimes curr_t]; %record the time, when exocytic event has happened. %Masha debug
            
   end

%% an endocytotic event may happen

    dice = random('unif',0,1); % dice=random number between 0 and 1
    if ( dice < dt/periodOut && curr_t ) % endocytosis will happen if dice<dt/periodOut (=0.05/0.6=0.08: 8% chance)  
       
      % site of endocytosis
      dice = random('unif',0,1); 
	  if ( dice < endo_m )                          % endo within window 
        dice = random('unif',0,1,N,N).*ves_prob;    % position in window
        [indx,indy] = find(dice==max(max(dice)));   % endo bin
      else                                          % endo outside window
        dice = random('unif',0,1,N,N).*ves_outprob; % position outside window
        [indx,indy] = find(dice==max(max(dice)));   % endo bin
      end
      % in case 2 bins get same random number
      if (length(indx)>1) 
        indx = indx(1); indy = indy(1);
      end


      % put this bin onto array with bins destined to endocytose
      % set negative endocytosis time for bins not yet full
      endo_cell.endo_time = [endo_cell.endo_time,-1]; 
      endo_cell.pick_time = [endo_cell.pick_time,curr_t];

      % adding proteins to vesicle
      endo_cell.vSNARE  = [endo_cell.vSNARE,   vSNARE(indx,indy)]; 
      endo_cell.Rec     = [endo_cell.Rec,      Rec(indx,indy)]; %Masha dev oct2012
      endo_cell.RecA    = [endo_cell.RecA,     RecA(indx,indy)]; %Masha dev oct2012
      endo_cell.Cdc42T  = [endo_cell.Cdc42T,   Cdc42T(indx,indy)];
      endo_cell.Cdc42D  = [endo_cell.Cdc42D,   Cdc42D(indx,indy)];
      endo_cell.BemGEF42= [endo_cell.BemGEF42, BemGEF42(indx,indy)];
      endo_cell.GDI42   = [endo_cell.GDI42,    GDI42(indx,indy)];
      endo_cell.BemGEF  = [endo_cell.BemGEF,   BemGEF(indx,indy)];      
      
      % recording the position of vesicle
      endo_cell.posx = [endo_cell.posx,indx*dx];
      endo_cell.posy = [endo_cell.posy,indy*dx];
      % recording the dx value at time vesicle was created
      endo_cell.dx   = [endo_cell.dx,dx];
      endo_cell.num  = endo_cell.num+1;
      

      % Radially 'squeeze' [Cdc42] distribution, set up variables for interpolation
      % 
      
      r_endo  = nVOutlen*dx/(pi^0.5);  % radius of area removed - old membrane size
                                       % area removed is (nVOutlen*dx)^2
      C = N*dx;
      hole_radial_endo_area_04;                         % call script 'hole' to compute new grid
                                   % in this case the new grid has a point
                                   % missing at the endocytosis point                                   

      % transform all concentrations so vesicle event is in the centre
      vesicle_position_center;
      
      % endo_squeeze - interpolates on the new grid
      %              - removes mass at vesicle from the membrane
      %                (scales it, to account for fewer grid points)
      %              - performs anotherscaling to keep within old max & min
      vSNARE   = endo_squeeze_radial ( vSNARE,   Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      Rec      = endo_squeeze_radial ( Rec,      Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2)); %Masha dev oct2012
      RecA     = endo_squeeze_radial ( RecA,     Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2)); %Masha dev oct2012
      Cdc42T   = endo_squeeze_radial ( Cdc42T,   Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      BemGEF42 = endo_squeeze_radial ( BemGEF42, Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      BemGEF   = endo_squeeze_radial ( BemGEF,   Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      Cdc42D   = endo_squeeze_radial ( Cdc42D,   Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      GDI42    = endo_squeeze_radial ( GDI42,    Fwarp,  x0, y0, x1, y1, N/2, N/2, N^2/(N^2-nVOutlen^2));
      clear xold yold
      
      % transform all concentrations back
      vesicle_position_restore;
      
      %% move endocytic patches after endocytosis

      dx_new  = dx*sqrt(N^2-nVOutlen^2)/N;
      C = N*dx_new;                    % move patches on new membrane size
      endo_cell = moving_patches_endo_radial(r_endo,indx*dx_new,indy*dx_new,endo_cell,C);  
      
      %%  

      % adjust size of membrane, and diffusion operator
      dx  = dx*sqrt(N^2-nVOutlen^2)/N;
      eta = ((Rnew_mult^3 * dx^3) - (((Rnew_mult * dx)-mem_depth)^3))/cyt_mult;
      Dxx = Dxx*(N^2-nVOutlen^2)/N^2;
      Hop = speye(N^2) - kron(speye(N),Dxx) - kron(Dxx,speye(N));
      
   
    end   

    % look for bins that are full, then have them wait 15 s
    ind = find( ...
        (endo_cell.vSNARE > patch_fill | curr_t-endo_cell.pick_time > patch_fill_give_up) ...
        & endo_cell.endo_time <= 0);
    endo_cell.endo_time(ind) = curr_t + endo_red_time;

    % look for bins that need to endocytose now
    ind = find((endo_cell.vSNARE>patch_fill | curr_t-endo_cell.pick_time ...
                  > patch_fill_give_up) & endo_cell.endo_time<=curr_t+0.5*dt & endo_cell.endo_time>=0);
    % if last condition false, that means bins need to turn red and wait 15 s
    
    if (~isempty(ind))  % a cell need to endocytose!
        % length(ind)
      endotimes=[endotimes curr_t];
        
      % only adding one vesicle to inner membrane despite the face that
      % length(ind) could be greater than one
      mydx2 = sqrt( dx2^2 + sum((nVOutlen*endo_cell.dx(ind)).^2 ) /N^2 );

%      fprintf(rec_count_time,'%5.4f %5.4f %5.4f %5.4f %5.4f\n',mean(mean(Rec)), sum(endo_cell.Rec(ind)), Recic, mean(mean(RecA)),sum(endo_cell.RecA(ind)) );
     
 
      % dump sink content to internal compartment
      vSNAREic = ((vSNAREic*N^2*dx2^2) + ...
          (sum(endo_cell.vSNARE(ind).* (endo_cell.dx(ind).^2) )) )/mydx2^2/N^2;  
      
      GDI42ic  = ((GDI42ic*N^2*dx2^2) + ...
           sum(endo_cell.GDI42(ind) .*((endo_cell.dx(ind).^2)) )  )/mydx2^2/N^2;     
      
      Cdc42Dic = ((Cdc42Dic*N^2)*dx2^2 + ...
            sum( endo_cell.Cdc42D(ind)                            .*((endo_cell.dx(ind).^2)) ) ...
          + sum( (endo_cell.Cdc42T(ind) + endo_cell.BemGEF42(ind)).*((endo_cell.dx(ind).^2)) ) ...
          )/mydx2^2/N^2;      
       %recemptor complex Rec sinking. If it needs to disappear, just comment
      %this line. The endosytic event will not change the receptor
      %concentration. Or can store in a different container, just to keep
      %track how much were internalized.
      Recic= ((Recic*N^2*dx2^2) + ...
          (sum(endo_cell.Rec(ind).* (endo_cell.dx(ind).^2) ))...
          + (sum(endo_cell.RecA(ind).* (endo_cell.dx(ind).^2) )) )/mydx2^2/N^2; %Masha dev dec2012
      
      % adding total Bem1 to cytoplasm      
      BemGEFc = BemGEFc + ...
          (eta*sum( (endo_cell.BemGEF42(ind) + endo_cell.BemGEF(ind)).*(endo_cell.dx(ind).^2)/(dx^2) )/(N^2)); 
      
      dx2 = mydx2;      
      
      % record lifetime, [vSNARE] & internalisation time of these bins
      record_lifetime = [record_lifetime (curr_t - endo_cell.pick_time(ind))];
      record_vSNARE   = [record_vSNARE   endo_cell.vSNARE(ind)];
      record_Rec = [record_Rec endo_cell.Rec(ind)];
      record_RecA     = [record_RecA     endo_cell.RecA(ind)];
      record_endot    = [record_endot    ones(1,length(ind))*curr_t];
      
      % take this bin off the array

      newind = setdiff([1:endo_cell.num],ind); % values in '[1:endo_cell.num]' = 1,2,3,...,endo_cell.num
                                               % and not in 'ind'
      endo_cell.endo_time = endo_cell.endo_time(newind);
      endo_cell.pick_time = endo_cell.pick_time(newind);
      
      endo_cell.vSNARE    = endo_cell.vSNARE(newind);
      endo_cell.Rec  = endo_cell.Rec(newind); %Masha dev oct2012
      endo_cell.RecA      = endo_cell.RecA(newind);     %Masha dev oct2012
      endo_cell.Cdc42T    = endo_cell.Cdc42T(newind);
      endo_cell.Cdc42D    = endo_cell.Cdc42D(newind);
      endo_cell.BemGEF42  = endo_cell.BemGEF42(newind);
      endo_cell.GDI42     = endo_cell.GDI42(newind);
      endo_cell.BemGEF    = endo_cell.BemGEF(newind);
      
      endo_cell.posx      = endo_cell.posx(newind);
      endo_cell.posy      = endo_cell.posy(newind);
      endo_cell.dx        = endo_cell.dx(newind);
      endo_cell.num       = endo_cell.num - length(ind);
 
    end
    
